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FOREWORD 

This report summarizes work performed during the first 
ten months of Contract NAS8 -27012, "Study of Propellant Dy- 
namics in a Shuttle -Type Launch Vehicle. " The work was con- 
ducted by Lockheed Missiles & Space Company, Inc., Huntsville 
Research & Engineering Center, for George C. Marshall Space 
Flight Center of the National Aeronautics and Space Administra- 
tion. The contract is administered under the direction of the 
Aero-Astrodynamics Laboratory, NASA-MSFC, with Mr. Larry 
Kiefling as Contracting Officer Representative. 
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SUMMARY 

This report describes a method and an associated digital computer program 
for evaluating the vibrational characteristics of large liquid -filled rigid wall 
tanks of general shape. A solution procedure was developed in which slosh 
modes and frequencies are computed for systems mathematically modeled as 
assemblages of liquid finite elements. To retain sparsity in the assembled 
system mass and stiffness matrices, a compressible liquid element formula- 
tion was incorporated in the program. 

The approach taken in the liquid finite element formulation is compatible 
with triangular and quadrilateral structural finite elements so that the analysis 
of liquid motion can be coupled with flexible tank wall motion at some future 
time. The liquid element repertoire developed during the course of this study 
consists of a two-dimensional triangular element and a three-dimensional 
tetrahedral element. 

A number of example problems were analyzed with the propellant dy- 
namics computer program. The results of these analyses compared favorably 
with known closed form solutions. Solutions computed for several finite element 
models of liquid sloshing in a rigid rectangular container are presented. 
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Section 1 
INTRODUCTION 

During this study, a method and an associated digital computer program 
were developed for computing the free undamped vibrational characteristics 
of large liquid -filled tanks of general shape. Solution procedures were formu- 
lated in which systems are mathematically modeled as assemblages of liquid 
finite elements. 

Unlike previous investigations (Refs. 1, 2, and 3), the methods developed 
during this study satisfy two fundamental requirements: (1) the procedure is 
suitable for analyzing large systems such as Space Shuttle vehicles, which 
typically involve thousands of degrees of freedom; and (2) the approach to liquid 
representation is compatible with triangular and quadrilateral structural finite 
elements, so that the analysis of liquid motion can be coupled with flexible tank 
wall motion at some future time. 

Solutions to eigenvalue problems resulting from large finite element 
networks can be obtained economically only if the assembled system mass and 
stiffness matrices exhibit sparse and/or banded characteristics. Most past 
investigations of finite element liquid sloshing problems that are reported in 
the literature are based on the assumption that fluid is incompressible. This 
assumption invariably yields system matrices that are completely full; conse- 
quently, the resulting computer programs are restricted to problems of moderate 
size (usually 200 degrees of freedom or less). In this study the effects of liquid 
compressibility are included in the finite element formulation, and the resulting 
system matrices exhibit the same sparse characteristics as an analogous three- 
dimensional solid finite element formulation. As a result, sparse and/or 
bandmatrix solution procedures (such as the matrix iteration method in- 
corporated in the Lockheed -Huntsville SNAP/Dynamics program, Ref. 4 can 
be used to solve the 0) MX - KX = 0 eigenvalue problem. 
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The mathematical model of a fluid system consists of liquid finite element 
representations that can easily be merged with solid shell finite element repre- 
sentations of tank walls. The element repertoire developed during the course 
of the study consists of a two-dimensional triangular element and a three- 
dimensional tetrahedral element. Particular care was taken to ensure complete 
compatibility between the liquid element formulations and the solid plate and 
membrane formulations incorporated in the SNAP/Dynamics program. As in 
SNAP, individual element energy matrices are formulated relative to imbedded 
intrinsic reference frames. This apparatus makes it easy to add new element 
formulations to the program, since it is necessary only to construct subrou- 
tines for computing the corresponding intrinsic mass and stiffness matrices, 
with element dimensions and fluid properties supplied through the calling 
sequence from general routines. 

The generalized coordinates used to characterize individual element 
energies are three -component displacement vectors associated with each 
element node. System mass and stiffness matrices are constructed by ex- 
pressing total system energies as the sum of individual element energies and 
invoking displacement compatibility among connecting elements through their 
common nodes. 
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Section 2 

FORMULATION OF THE EIGENVALUE PROBLEM 

2.1 SYSTEM EQUATIONS OF MOTION 

The kinetic and potential energies of a system may be expressed in 
matrix form as: 


T = J <f>* M 4> 

and ( 1 ) 

1 * 

V K $ , 

where $ is the generalized coordinate vector, and M and K are the mass 
and stiffness matrices of the system. In the absence of dissipative effects and 
external forces, the Lagrange equations yield 


M <3> + K = 0 . 


(2) 


Assuming solutions are in the form of <J> = X cosCOt, Eq. (2) is reduced to the 
usual linear vibrational eigenproblem 


00 2 MX-KX=O . 


(3) 
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In this study, node point displacement components are used explicitly 
as system generalized coordinates, so that for a network composed of a total 
of n system nodes 



(4) 


where v 1 is a vector containing the three motion components of node i. 

The total system kinetic and potential energies may be expressed as 


T 


N 

Ex.. 

j=i j 


and 


V 



» 


(5) 




where and represent the kinetic and potential energies of the j liquid 
element in terms of the system generalized coordinates, and N is the total 
number of finite elements in the network. 


The kinetic and potential energies of liquid element j are expressed as 


1 • * • 

T. = T <p. M. (P. t and 
3 2 3 3 3 


V. 

J 




* 


K. cp. 
3 3 


9 


( 6 ) 
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where M. and K. are the mass and stiffness matrices of the element in terms 
J J 

of element generalized coordinates. In the propellant dynamics computer pro- 
gram, element energies are expressed in terms of explicit nodal displacements; 
so that, where u 1 is a vector containing the three motion components of node i 
in directions parallel to the intrinsic reference frame axes of element j, 





(7) 


In Eq. (6), m is number of nodes interconnected by element j. 

To transform element energy matrices M. and K. into system coordinates, 

J J 

a coordinate transformation is performed such that 




R. <f> 
3 


( 8 ) 


Substituting Eq. (8) into Eqs. (6) yields 


where 


— 1 — ; 

T . = j- <p M. <P , and 
J 2 J 

V. =4$ K. $ f 

j 2 j 


M. = R. M. R., and 
J J J J 


K. = R * K. R. . 
J J J J 


(9) 


( 10 ) 
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Substitution of Eqs. (9) into Eqs. (5) yields the total system kinetic and potential 
energies in terms of the desired coordinates. 

For purposes of specifying node restraints and system boundary condi- 
tions, unique reference frames are associated with each system node, and the 
directions of node displacements are defined relative to these nodal reference 
frames. Consequently, the transformation matrix IF includes the effects 
of: (1) the orientation of the intrinsic reference frame of element j; and (2) 
the orientations of the nodal reference frames associated with nodes 1 through 
m. 


The procedure described above is identical to that incorporated in SNAP/ 

Dynamics, and the sparsity of the resulting system matrices is the same as 

would be obtained in a solid finite element formulation. However, if the element 

matrices, M. and K., were formulated assuming liquid incompressibility, 

J 1 

modification of the system matrices produced by Eqs. (5) would be required 
before solution to the eigenproblem could be attempted. Expressions would 
be introduced to ensure that the entire liquid volume would undergo zero volume 
change. These expressions would necessarily relate the relative motions of 
all network nodes that appear on the surfaces of the finite element model. Conse- 
quently, sparsity in the assembled system matrices would be lost, and those 
rows and columns associated with surface nodes would be completely full. This 
procedure would drastically restrict the size capability of the resulting computer 
program. By incorporating the effects of compressibility in the element formu- 
lations, the system matrices produced by Eqs. (5) can be used directly, and 
sparse and/or bandmatrix procedures can be used to solve the resulting eigen- 
value problem. 


2.2 LIQUID FINITE ELEMENTS 

The liquid is assumed to be inviscid and compressible. The effect of 
small variation in mass density caused by the compressible assumption is 
neglected in formulating element kinetic energies. 


2-4 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC-HREC D225823 


The intrinsic reference frame associated with a tetrahedral liquid ele- 
ment is illustrated on Fig. 1. The origin of the reference frame coincides with 
node 1; intrinsic axis 1 extends from node 1 to node 2; axis 2 lies in the plane 
of nodes 1, 2, and 3 and is directed so that node 3 has a positive 2 coordinate. 
The relative displacements q^, q^, . . . , and q^ as shown on Fig. 1 completely 
define the deformation of the element. As the element moves with the liquid 
network, the intrinsic reference frame travels with node 1. The motion of 
the reference frame is represented by q^, ...» and q^, as illustrated on 

Fig. 1. Liquid element mass and stiffness matrices are derived in terms of 
the twelve coordinates q^ through q^ and q^ through q^. 

A linear deformation field is assumed for the element; it is expressed as 


w . 1 

1 

X, X-, 

x_ 


1 1 


1 , 2 

3 


W 2 [ 

1 

X 1 x 2 

x 3 

A Q, or 

W 3 ) 

1 

X 1 X 2 

X 3_ 



w = 

X A Q , 



( 11 ) 


( 12 ) 


where w defines the motion of a particle on the interior of the element with 
intrinsic coordinates (x^, x^, x^) and Q is a vector of the twelve coordinates 
q ^ through q^ and q^ through q^. The terms of the coefficient matrix A 
are determined according to the physical dimensions of the element. 

The kinetic energy of a liquid element is expressed as the following 
volume integral 


T 



volume 


(13) 
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Fig. 1 - Intrinsic Reference Frame and Generalized Coordinates 
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where p is the mass density of the liquid. Substituting Eq. (12) into Eq. (13) 
yields the element kinetic energy as a quadratic form in terms of the intrinsic 
generalized coordinates, so that 


T 



M Q 


9 


(14) 


where M represents the intrinsic element mass matrix. 

The total potential energy of a liquid element is expressed as the summa- 
tion of two energies: (1) gravitational potential; and (2) dilatational potential. 
The potential energy incurred by gravity considerations is expressed by the 
following surface integral 


? G = ’ < 15 > 

surface 

where g n is the component of the gravitational acceleration vector normal to 
the surface differential, dS, and £ is the component of surface motion normal 
to the surface differential. Dilatational energy is expressed as 


V D = \ «J e 2 dv , (16) 

volume 

where K is the liquid bulk modulus and 0 represents the volume change of the 
element. Substituting into Eqs. (15) and (16) the normal component of surface 
motion, £, and the element dilatation, 0, in terms of the assumed displacement 
field represented by Eqs. (12) yields the element potential energies as quadratic 
forms in terms of the intrinsic generalized coordinates, so that 
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V 


G 


1 * ~ 

2 Q Kq Q , and 






Q . 


The total element potential energy is 


(17) 


V = V G + V D * ° r 

V = | Q" K Q , 


(18) 


where K represents the intrinsic element stiffness matrix and is equal to the 

summation of K„ and K^. 

Lr JJ 


A transformation of coordinates is performed so that the element energy 
matrices appearing in Eqs. (14) and (18) are expressed in terms of explicit 
nodal displacements. The resulting matrices are then used to form the as- 
sembled system mass and stiffness matrices as described in the previous 
section. 


Energy expressions for the two-dimensional triangular element included 
in the program are derived in a manner similar to the tetrahedral element. 
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Section 3 

COMPUTER PROGRAM 

During this study a general purpose digital computer program was de- 
veloped for computing the slosh modes and frequencies of finite element repre- 
sentations of liquid filled rigid wall tanks of general shape. Some of the features 
of the program are as follows: 


• The element repertoire consists of a two-dimensional triangular 
element and a three-dimensional tetrahedral element. 

• All operations involved in forming element and system matrices 
are performed in double precision. 

• A moderate-size eigensolution routine is incorporated that 
will accommodate a? M X - KX = 0 eigenproblems of up to 
approximately 75 degrees of freedom in double precision and 
100 degrees of freedom in single precision. 

• System nodes may be resequenced, giving the user maximum 
control over the bandwidth of system matrices. 

• Individual reference frames may be assigned to the system 
nodes. Nodal motion components are given relative to these 
reference frames. 

• Automatic network generators are incorporated to facilitate 
the input of node reference frame assignments, node restraint 
conditions, node position coordinates, and element definitions. 

• SC 4020 plotting of undeformed and deformed element networks. 

• Intrinsic local coordinate systems are used for individual elements. 


Like SNAP/Dynamics, the program is arranged in a modular fashion to 
provide a convenient framework for modifying and expanding various parts of 
the program. New element formulations can be included in the element repertoire 
without altering any of the existing code. The module for computing eigensolu- 
tions can also be easily replaced. 
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The program is coded in Fortran V and designed for use on the NASA- 
MSFC Univac 1108 Exec VIII system. 
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Section 4 
RESULTS 

Solutions computed for liquid sloshing in a rigid rectangular container 
are presented. The solutions are compared with known closed-form solutions. 
The geometric dimensions of the liquid volume are 1 meter wide, 2 meters 
long and 1 meter deep. The mass density, bulk modulus and equivalent gravita- 

tionsl acceleration of the liquid are 1000 kg/m^, 2 x 10^ kg/m-sec^ and 9.807 

/ 2 

m/sec , respectively. 

Results for six finite element models of the same liquid volume are pre- 
sented. Three of the networks are two-dimensional and are composed of tri- 
angular elements. The two-dimensional models, which are referred to as 
2D-a, 2D-b, and 2D-c, are illustrated on Fig. 2. The other three models are 
three-dimensional and are made up of tetrahedral elements. The three-di- 
mensional models, which are referred to as 3D-a, 3D-b, and 3D-c, are shown 
on Fig. 3. The number of elements contained in each model and frequency com- 
parisons between closed-form solutions are given in Table 1. Plots illustrating 
the slosh modes computed with the two-dimensional models are given in Fig. 4. 
Plots illustrating the modes computed with the three-dimensional models are 
given in Fig. 5. 
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COMPARISON OF FREQUENCIES 
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3D-c 96 0.589449 0.815962 
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Mode 1, Frequency = 0.595833 cps 



Mode 2, Frequency = 0.848202 cps 

Fig. 4 - (Continued) — Case 2D-b 
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Mode 4, Frequency = 1.02322 cps 
Fig. 4 - (Continued) — Case 2D-b 
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Mode 1, Frequency = 0.598934 cps 



Mode 2, Frequency = 0.889266 cps 
Fig. 4 - (Continued) — Case 2D-c 
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Mode 3, Frequency = 1.05944 cps 



Mode 4, Frequency = 1.08128 cps 
Fig. 4 - (Continued) — Case 2D-c 
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Mode 1, Frequency = 0.588683 cps 



Mode 2, Frequency = 0.709116 cps 


Fig. 5 - Free Surface Mode Shapes — Case 3D-a 
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Mode 1, Frequency = 0.589449 cps 



Mode 2, Frequency = 0.815962 cps 


Fig. 5 - (Continued) — Case 3D-c 
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